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ABSTRACT 

We study the orbital evolution and accretion history of massive black hole (MBH) 
pairs in rotationally supported circumnuclear discs up to the point where MBHs form 
binary systems. Our simulations have high resolution in mass and space which, for 
the first time, makes it feasible to follow the orbital decay of a MBH either counter- 
or co-rotating with respect to the circumnuclear disc. We show that a moving MBH 
on an initially counter-rotating orbit experiences an "orbital angular momentum flip" 
due to the gas-dynamical friction, i.e., it starts to corotate with the disc before a 
MBH binary forms. We stress that this effect can only be captured in very high 
resolution simulations. Given the extremely large number of gas particles used, the 
dynamical range is sufficiently large to resolve the Bondi-Hoyle-Lyttleton radii of 
individual MBHs. As a consequence, we are able to link the accretion processes to the 
orbital evolution of the MBH pairs. We predict that the accretion rate is significantly 
suppressed and extremely variable when the MBH is moving on a retrograde orbit. It 
is only after the orbital angular momentum flip has taken place that the secondary 
rapidly "lights up" at which point both MBHs can accrete near the Eddington rate 
for a few Myr. The separation of the double nucleus is expected to be around < 10 
pc at this stage. We show that the accretion rate can be highly variable also when 
the MBH is co-rotating with the disc (albeit to a lesser extent) provided that its 
orbit is eccentric. Our results have significant consequences for the expected number 
of observable double AGNs at separations of < 100 pc. 

Key words: black hole physics - hydrodynamics - galaxies: starburst - galaxies: 
evolution - galaxies: nuclei 



1 INTRODUCTION 

In recent years, compelling observational evidence has in- 
dicated that massive black holes (MBHs) are ubiquitous in 
nuclei of local bright galaxies (Richstone et al. 1998; Mar- 
coni et al. 2004; Shankar et al. 2004; Ferrarese et al. 2006; 
Ferrarese 2006; Decarli et al. 2007). According to the struc- 
ture formation paradigm, galaxies often interact and merge 
as their dark matter halos assemble in a hierarchical fashion, 
and the MBHs, incorporated through mergers, are expected 
to grow, evolve and pair with other MBHs (Begelman Bland- 
ford & Rees 1980; Volonteri, Haardt, & Madau 2003). The 
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formation of MBH pairs thus appears to be an inevitable 
and natural consequence of galaxy assembly. 

A large number of merging systems such as the lumi- 
nous infrared galaxies (LIRGs) hosts a central rotationally 
supported massive (up to 10 10 M Q ) gaseous disc extending 
on scales of ~ 100 pc (Sanders & Mirabel 1996; Downes & 
Solomon 1998; Davies et al. 2004). These discs may be the 
end-product of gas-dynamical, gravitational torques excited 
during the merger, when a large amount of gas is driven into 
the core of the remnant (Mayer et al. 2007). Inside a massive 
self-gravitating disc, a putative MBH pair can continue its 
dynamical evolution and, possibly, accrete gas producing a 
double AGN (Kocsis et al. 2005, Dotti et al. 2006a). 

The possibility of double accretion processes during the 
different stages of a galaxy merger is still a matter of debate. 
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From an observational point of view, only few tens of dou- 
ble quasars are known (with separations < 100 kpc; see e.g. 
Foreman, Volonteri & Dotti 2008 and references therein). 
There are only a few well studied cases of ongoing merg- 
ers at sub-galactic scales in which each companion is an 
AGN. Examples of such systems include NGC 6240 (Ko- 
mossa et al. 2003), Arp 229 (Ballo et al. 2004), and Mrk 
463 (Bianchi et al. 2008). Other two dual AGN candidates 
with separation ~ 1 kpc have been identified spectroscopi- 
cally by Comerford et al. (2008). On parsec scale, only one 
resolved active MBH binary has been found in the nucleus 
of the elliptical galaxy 0402+369 (Rodriguez et al. 2006), 
and the existence of two more sub-parsec MBH binary can- 
didates has been suggested (OJ287, Valtonen et al. 2007, 
and SDSSJ092712.65+294344.0, Bogdanovic, Eracleous & 
Sigurdsson 2008; Dotti et al. 2008). 

Accretion events during the last phases of a MBH bi- 
nary (for sub-parsec separations) have been studied in de- 
tail by several authors (Armitage & Natarajan 2002, 2005; 
Hayasaki, Mineshige, & Sudou 2007; Hayasaki, Mineshige, & 
Ho 2008; Cuadra et al. 2009). On larger scales, the possibility 
of gas accretion during the MBH pairing has been recently 
studied using numerical simulations in the context of galaxy 
mergers (e.g., Kazantzidis et al. 2005; Springelet al. 2005; Di 
Matteo et al. 2005; Hopkins et al. 2005; Hopkins et al. 2006). 
In these works, the authors study a galaxy-galaxy collision 
on spatial scales of ~ 100 kpc, and the number of particles 
used to model the two galaxies does not allow them to re- 
solve the influence radii of the two pairing MBHs. This lack 
of resolution prevents the authors from accurately predict- 
ing the accretion rates onto the two MBHs. Given the limits 
of the currently available supercomputers, such simulations 
would be prohibitively time consuming. Therefore, no sim- 
ulations resolving the Bondi-Hoyle-Lyttleton radii during a 
complete galaxy-galaxy merger have been published to date. 

We have run a suite of high resolution N- 
body/hydrodynamical simulations of MBH pairing in cir- 
cumnuclear discs. In this paper we discuss the dynamical 
evolution of the MBH pairs up to a point where they form 
a binary (i.e., up to a point where the mass enclosed inside 
the MBH orbit is smaller than the sum of the MBH masses) . 
We follow the MBH dynamics from an initial separation of 
50 pc down to sub-parsec separation while resolving -Rbhl 
of the two pairing MBHs. We predict the accretion rate onto 
the two MBHs as a function of the geometry of the merger. 
As we demonstrate in this paper, this approach allows us 
to study the possibility that a galaxy merger event may be 
associated with a single or double active galactic nucleus 
depending on the configuration of the initial MBH orbits. 
It also allows us to constrain the fate of the MBH binary 
(i.e., coalescence vs. stalling). The evolution of the MBH 
spins and the fate of the newly formed MBH binary will be 
discussed in two forthcoming papers. 



2 SIMULATION SETUP 

We follow the dynamics of MBH pairs in nuclear discs us- 
ing numerical simulations run with the N-Body/SPH code 
GADGET (Springel, Yoshida & White 2001), upgraded to 
include the accretion physics. The main input parameters of 
our simulations are summarized in Table 1. 



Table 1. Run parameters 



run 


prograde ? 


e; 


7 


HPC 


yes 







HPE 


yes 


0.7 


5/3, "hot" 


HRE 


no 


0.7 




CPC 


yes 







CPE 


yes 


0.7 


7/5, "cold" 


CRE 


no 


0.7 





The initial conditions of our runs are set following the 
procedure discussed in Dotti, Colpi & Haardt (2006b) and 
Dotti et al. (2007). In our models, two MBHs are placed 
in the plane of a gaseous disc, embedded in a larger stellar 
spheroid. The gaseous disc is modeled with 2,353,310 parti- 
cles, has a total mass JVfDisc = 10 8 Mq, and follows a Mestel 
surface density profile 

M disc 



E(7?) = 



(1) 



27ri?discfl 

where R is the radial distance projected into the disc plane 
and -Rdisc = 100 pc is the size of the disc. The disc is rota- 
tionally supported in R and has a vertical thickness of 10 
pc. Initially, the gaseous particles are distributed uniformly 
along the vertical axis. The internal energy per unit mass of 
the SPH particles scales as: 



u(R) = KR- 2/:i , 



(2) 



where K is & constant denned so that the Toomre parameter 
of the disc, 

KC S 



Q = 



ttGE' 



(3) 



is > 3 everywhere, preventing the fragmentation of the disc 
(the average value of Q over the disc surface is « 10). In 
equation (3) k is the local epicyclic frequency, and c s is the 
local sound speed of the gas. Gas is evolved adiabatically as- 
suming a polytropic index 7 = 5/3 or 7 = 7/5. In the former 
case, the runs are denoted by "1" and are termed "hot" as 
the temperature is proportional to a higher power of density 
than in the latter class of runs ("cold" cases, runs denoted 
by "2"). The hot case corresponds to adiabatic monoatomic 
gas, while the cold has been shown to provide a good ap- 
proximation to a solar metallicity gas heated by a starburst 
(Spaans & Silk 2000; Klessen, Spaans, & Jappsen 2007). 

The spheroidal component (bulge) is modeled with 10 5 
collisionless particles, initially distributed as a Plummer 
sphere with mass density profile 

M / 2 ^ ~ 5/2 

JWBulgo 



p(r) 



47T b' A 



1 + ¥ 



(4) 



where b (= 50 pc ) is the core radius, r the radial coordinate, 
and MBuig C (= 6.98MDi sc ) is the total mass of the spheroid. 
With such a choice, the mass of the bulge within 100 pc is 
five times the mass of the disc, as suggested by Downes & 
Solomon (1998). 

The two MBHs (Mi and M 2 ) are equal in mass (M B h = 
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4xl0 6 M ). The initial separation of the MBHs is 50 pc. Mi, 
called primary for reference, is placed at rest at the centre 
of the circumnuclear disc. The secondary (denoted as M2) is 
either corotating with the gaseous disc on a circular orbit, or 
moving on an eccentric counterrotating or corotating orbit 
with respect to the circumnuclear disc. For eccentric orbits, 
the initial ratio between the absolute values of the radial 
and tangential velocity of M2 is « r ad/«tan = 3. This ratio 
sets the eccentricity of the first orbit of M2 toe~ 0.7. Given 
the large masses of the disc and the bulge, the dynamics of 
the moving MBH (M2) is uneffected by the presence of Mi 
until the MBHs form a gravitationally bound system. The 
subsequent evolution of the MBH binary will be discussed 
in a forthcoming paper. 

We evolve our initial composite model (bulge, disc and 
Mi) for « 3 Myrs, until the bulge and the disc reach equi- 
librium, as described in Dotti et al. (2006b; 2007). The disc, 
having initially finite thickness and homogeneous vertical 
density distribution, is allowed to relax to an equilibrium 
configuration along the z-axis. The new distribution of the 
gas in the disc in the z direction becomes non-uniform once 
the equilibrium configuration is reached. The new disc has 
a vertical thickness of ~ 8 pc (defined as containing 90% 
of the mass in the vertical direction at every radius) This 
thickness is independent of the radius and is « 4/5 of the 
initial un-relaxed value. 

We allow the gas particles to be accreted onto the MBHs 
if the following two criteria are fulfilled: 

• the sum of the kinetic and internal energy of the gas 
particle is lower than a times the absolute value of its grav- 
itational energy (all the energies are computed with respect 
to each MBH) 

• the total mass accreted per unit time onto a MBH ev- 
ery timestep is lower than the accretion rate corresponding 
to the Eddington luminosity (L^dd) computed assuming a 
radiative efficiency (e) of 10%. 

The parameter a is a constant defining how much a 
particle has to be bound to a MBH to be accreted. We set 
a = 0.3. Note that due to the nature of the above crite- 
ria, the gas particles can accrete onto the MBHs only if the 
Bondi-Hoyle-Lyttleton radius, defined as 



is resolved in the simulations. In Eq. (5), M is the mass 
of the MBH, v TO i is the relative velocity between the gas 
and the moving MBH, and c s is the sound speed of the gas. 
M is evolved by adding at each timestep 1 — e times the 
total mass of the accreted gas particles. If the amount of 
bound gas particles implies super-Eddington accretion, the 
code randomly selects a subsample of the particles within 
7?bhl to be accreted, so as to prevent accretion from reach- 
ing super-Eddington rates. This simplified treatment of the 
black hole feedback does not capture effects such as the ac- 
cretion driven wind, and its effect on the environment. Be- 
cause the MBHs accrete only a small fraction of their ini- 
tial masses (< 10%) during our simulations, AGN feedback 
is not expected to modify the global properties of the cir- 
cumnuclear disc. Nevertheless, such feedback can remove a 
fraction of the gas a few parsecs away from MBHs accreting 



near the Eddington limit. This could decrease the accretion 
rate onto the MBHs and the efficiency of MBH pairing after 
they form a binary system. The accurate implementation of 
a self-consistent radiative feedback is beyond the scope of 
this work, and will be the subject of our future studies. 

The number of neighboring SPH particles adopted in 
our simulation is iV no igh = 32. The resulting spatial reso- 
lution of the hydrodynamical force in the highest density 
regions is « 0.1 pc. In order to prevent numerical errors due 
to the density-dependent effective resolution, we set both 
the gravitational softening of the gaseous particles and that 
of the MBHs to 0.1 pc. With this spatial resolution we can 
resolve the vertical scale of the disc and, as a consequence, 
we can calculate angular momentum transport due to disc 
self-gravity (see, e.g., Lodato & Rice 2004; Nelson 2006). 
We can also resolve the influence radius of Mi (« 1 pc) and 
M2, a condition necessary to study the gas accretion onto 
MBHs t. 

In addition, our accretion scheme is such that linear mo- 
mentum is conserved. This includes automatically the brak- 
ing exerted by the accretion flow onto the moving MBH 
(Fbhl = Mv re i, where M is the accretion rate). The only 
viscosity term in the simulations is the shear reduced version 
(Balsara 1995; Steinmetz 1996) of the standard Monaghan 
and Gingold (1983) artificial viscosity. 



3 RESULTS 

3.1 Orbital evolution 

In this Section we discuss the orbital evolution of the pairs 
in six different runs, until the MBHs form a binary. The re- 
lationship between the orbital properties of the MBHs and 
their accretion properties is discussed in Section 3.2. The 
run names and parameters are listed in Table 1. The letters 
used to name the runs encode the type of the run with the 
first, second and third letter corresponding to hot vs. cold, 
prograde vs. retrograde and eccentric vs. circular cases, re- 
spectively. For example, CRE corresponds to a secondary 
MBH on an eccentric retrograde orbit in a cold disk. 

3.1.1 Prograde orbits 

The pairing process of MBHs on circular orbits (runs HPC 
and CPC) brings the two MBHs down to separations < 5 
pc in less than 10 Myrs at which point they form a binary. 
The MBHs orbital decay proceeds without any observable 
eccentricity growth. Because of the low relative velocity be- 
tween the MBH and the gas, M2 interacts efficiently with 
the gas out to a distance of a few parsecs. In this case, the 
dynamical evolution of the pair is almost indistinguishable 
from lower (parsec scale) resolution simulations (Escala et 
al. 2005; Dotti et al. 2006b; Dotti et al. 2007), before the 
MBHs bind in a binary. 

The upper panel of Fig. 1 shows the separation between 
the two MBHs as a function of time for M2 moving initially 
on a prograde eccentric orbit inside the hot disc (7 = 5/3, 

t Note that i?BHL of A^2 depends on the phase of its orbit. That 
is, i?BHL depends on the relative MBH-gas velocity, which is 
different in the corotating and counterrotating cases. 
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run HPE). The two MBHs reach a separation of the order 
of < 10 pc in less than 5 Myr. The MBH pair loses memory 
of its initial eccentricity in the early phases of the orbital 
decay. Such circularization is due to the variation in direc- 
tion and in magnitude of dynamical friction acting on M2 
during an orbital period. Dynamical friction tends to slow 
down the secondary near the periastron and to accelerate it 
at the apoastron (see Dotti et al. 2006b; Dotti et al. 2007 
for a detailed discussion). The orbital decay and circular- 
ization observed in these high resolution runs are slightly 
faster than in the low resolution simulations presented in 
Dotti et al. (2006b) and Dotti et al. (2007). The increased 
resolution allows to better resolve the interaction between 
the gas environment and M2 when it is still moving on ec- 
centric orbits and the relative velocity between M2 and the 
gas particles is high. This effect of the increased resolution 
is more prominent for retrograde orbits, where the relative 
velocities are higher. We postpone a detailed discussion of 
this effect to Section 3.1.2. 

The thermodynamics of the gas also has impact on the 
circularization process. The dynamical interaction of the 
moving MBH of mass M2 with the gas is stronger in a cold 
disc. This is especially so when the secondary is corotating 
with the disc. The reason for this is that the relative veloc- 
ity between the secondary and the disc material is small in 
this case (at least for circular orbits) and the Bondi-Hoyle- 
Lyttleton radius reduces to the Bondi radius which is pro- 
portional to cj 3 , where c s is the sound speed. It is for this 
reason that the orbital decay and the circularization in a 
cold disc (run CPE) proceeds faster than in the hot case. 
In our simulations dynamical friction is efficient in reducing 
the eccentricity down to values < 0.1, for both cold and hot 
cases. 

3.1.2 Retrograde orbits 

For the retrograde cases (runs HRE and CRE) we find that 
the two MBHs also reach a separation of < 10 pc in less than 
5 Myr, as illustrated in the upper panels of Fig. 2 and Fig. 3. 
The middle panels of Fig. 2 and Fig. 3 show the evolution 
of the z-component of the orbital angular momentum L z of 
M2 normalized to its initial value. The (initially negative) 
angular momentum grows very fast during the first Myr. 
As long as M2 is counter-rotating, the MBH-disc interac- 
tion brakes the secondary even at the apocentre, because 
the MBH velocity and disc flow are anti-aligned there. The 
increase of the orbital angular momentum of M2 is further 
facilitated by the fact that, while the orbit decays, the MBH 
interacts with progressively denser regions of the disc closer 
to the primary. Eventually, at about 2 Myr, the sign of the 
orbital momentum of the secondary changes. The dynamical 
friction process is the ultimate cause of this orbital angular 
momentum flip. After the orbital angular momentum flip 
has taken place and the secondary enters a co-rotating orbit 
with respect to the disc, the angular momentum continues 
to grow up to L/Lq ~ 0.1. 

The orbital decay of a retrograde MBH (runs HRE and 
CRE), and the orbital angular momentum flip process are 
dependent on the numerical resolution. We stress that the 
orbital angular momentum flip is not seen in lower resolution 
simulations. This resolution dependence 

has a simple dynamical explanation. The perturbation 
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2 4 

t [Myr] 

Figure 1. Run HPE. Upper panel: MBH separation as a func- 
tion of time. Middle panel: time evolution of the orbital angular 
momentum of M2 (L) normalized to its initial value (Lq). The 
thin-dashed horizontal line marks L = 0. Lower panel: Eddington 
accretion ratio as a function of time. Red and blue lines refer to 
M\ and M2, respectively. 



of the orbit of a single gas particle that is caused by M2 
depends on the relative velocity between the gas and the 
MBH and their relative separation. Only those gas particles 
that get perturbed can contribute to the orbital decay of 
the MBH and lead to its orbital angular momentum flip. A 
natural scale lenght for the gravitational interaction between 
the MBH and gas particles is -Rbhl- 

Gas particles can efficiently exchange angular momen- 
tum with the hole only if they pass closer to the MBH than 
a few -Rbhl- 

As long as the MBH is counter-rotating, its relative 
velocity with respect to the disc is v Ye \ > 200 km s _1 > 
c s , which corresponds to J?bhl J; 1 pc. Therefore, higher 
numerical resolution is required to model counter-rotating 
cases than the co-rotating ones. Given the high resolution 
of our simulations, the gravitational interaction is accurately 
computed down to scales of ft; 0.1 pc, allowing us to correctly 
model the dynamical evolution of counter-rotating MBHs. 

3.2 Accretion history 

In run HPC the Eddington ratio /Edd = M /M^dd of the pri- 
mary is « 1 at the beginning of the simulation and decreases 
monotonically during the pairing process down to ~ 0.7 at 5 
Myr. The Eddington ratio of the primary decreases steadily 
because the orbital motion of the secondary heats up the 
gas in the vicinity of the primary. The increase in the gas 
temperature decreases the Bondi radius of the primary and 
evacuates gas from the central regions of the disc. Both of 
these effects reduce the accretion rate onto the primary. This 
effect is present only when the polytropic index 7 = 5/3, as 
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Figure 2. Same as Fig. 1 for run HRE. In the middle panel 
the initial angular momentum is considered negative since M2 is 
initially orbiting on a retrograde orbit. 




Figures 1, 2 and 3 allow to compare directly the dy- 
namical properties of the two MBHs to the accretion rates 
in runs HPE, HRE, and CRE, respectively. Top panels show 
the evolution of the MBH separation. Middle panels present 
the time dependence of the orbital angular momentum of 
Mi. The bottom panels show /Edd for the primary (red lines) 
and secondary (blue lines) MBHs. 

In runs HPE and HRE, the primary accretes at an av- 
erage /Edd ~ 0.7, slightly decreasing with time. The ac- 
cretion rate of the secondary evolves differently. Its accre- 
tion behaviour depends on whether it is on a co-rotating or 
counter-rotating orbit. In the co-rotating case (run HPE), 
the average /Edd of the secondary is w 0.4 during the first 5 
Myr and its accretion history can be divided into two phases: 

i) fo r t ^ 2.5 Myr the circularization process is efficient 
while /Edd ~ 0.3 on average and shows strong variability 

ii) for t > 2.5 Myr, M2 moves on a quasi-circular orbit, 
and the relative velocity between M2 and the gaseous disc 
is reduced. In this phase, /Edd ~ 0.45. 

In the counter-rotating case (run HRE), we can still 
distinguish two phases albeit with some important differ- 
ences: 

i) for t < 3 Myr, M2 is counter-rotating (L z < 0), and 
/Edd ~ 0.05 -r 0.1. This level of accretion is lower than in 
the corresponding first orbital stage in the co-rotating case 
described above. Moreover, the variability in the accretion 
rate onto the secondary MBH is now significantly larger. 

ii) At t ~ 3 Myr the orbital angular momentum flip 
takes place. Thereafter, AI2 accretes at /Edd ~ 0.7. In this 
second stage, the main difference between the co-rotating 
and the counter-rotating cases is that, in the latter, M2 ex- 
hibits a rapid transition/increase in the accretion rate (cf. 
bottom panels in Fig. 2 (counter-rotating case) and Fig 1 
(co- rotating case). In the next subsection we present the 
Fourier analysis of the accretion history of M%, in order 
to assess whether the variability of /Edd is physical or it 
depends on numerical fluctuations due to finite number of 
particles used. 

In run CRE, M2 has a similar accretion history to the 
one observed in run HRE, but the timescale for the orbital 
angular momentum flip is shorter. Lower 7 corresponds to 
colder and denser gas and, thus, more efficient dynamical 
friction. This makes M2 corotate in less than 2 Myr. The 
primary has an average /Edd ~ 0.9, which is higher than in 
runs HPE and HRE. The physical reason for the enhance- 
ment in the accretion rate in the cold case is provided in the 
discussion. 



Figure 3. Same as Fig. 1 and Fig. 2 for run CRE. 



in that case the increase of the gas temperature with density 
is larger than in the cold case where 7 = 7/5 (run CPC). 
In the cold run, the lower value of 7 allows for the pairing 
of the MBHs while preventing a strong heating of the gas. 
In both runs (HPC and CPC), M2 has an accretion history 
similar to Mi. That is, MBHs on circular corotating orbits 
have low relative velocities with respect to the gas which 
leads to high and steady accretion rates. 



3.2.1 Fourier analysis of accretion fluctuations in run 
HRE 

As discussed above, strong fluctuations of /Edd are a distinc- 
tive feature of counter-rotating (and eccentric co-rotating) 
MBHs. The variability is higher for counter-rotating MBHs 
and /Edd < 0.1 (see Fig. 2 and 3). These values of /Edd cor- 
respond to less than a few hundred SPH particles accreted 
every 2.5 x 10 4 yr which is the time resolution in Fig. 2 and 
3. Such a small number of particles is < 10 x Af no i g h, and 
can drop down to when the -Rbhl of M2 is not resolved. 
As a consequence, values of /Edd < 0.1 should be considered 
order of magnitude estimates. 
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frequency 

Figure 4. Power spectrum P of /Edd fluctuations of M2 in run 
HRE, until the MBH is counter— rotating with respect to the cir- 
cumnuclear disc. Solid black line: P of the original /Edd- Dot- 
dashed blue line: P of the larger scale/noise free estimate of the 
/Edd discussed in the text. P has been normalized to its value 
at v = 0, v is in units of the average between v\ = 1/ri and 
V2 = 1/t2. Short dashed green line marks v = u\, long dashed 
red lines shows v = V2 and its higher harmonics. The dashed re- 
gions show the values of v compatible with the peaks of P (i.e. 
the frequencies of the peaks ±1/2 !^ m i n ). 



In order to check whether the variability of /Edd is due 
to finite sampling of the medium or if it has a physical mean- 
ing, we performed the Fourier analysis for run HRE*. We 
emphasise that run HRE is the run with the lowest /Edd, 
i.e., most affected by the numerical errors. Therefore our 
conclusions will be conservative and robust. 

The black solid line in Fig. 4 refers to the power spec- 
trum (P) of the fluctuations of /Edd of M2, for run HRE. 
The frequency (y) unit is denned as the average between 
v\ = 1/n and V2 = I/V2, where t\ and T2 are the periods 
of the first and second orbit of M%. P has been normal- 
ized to P[y = 0). The frequency resolution of our results is 
Vmin = 1/ At, where At « 3 Myr corresponds to the amount 
of time during which M2 orbit is retrograde. § 

Fig. 4 proves that P is not compatible with white noise 
as expected for fluctuations due only to numerical sampling. 
Moreover, we can distinguish four peaks in our distribution. 
The first peak at v = is a simple consequence of the phys- 
ical condition /Edd > 0. The remaining three peaks contain 
more physically useful information. Given our frequency res- 
olution, the frequency of the second peak is compatible with 
v\, while the third and fourth peaks are at the frequencies 

■f In run CRE M2 is counter-rotating for less than 2 orbits, so 
the Fourier analysis in this case would not be meaningful. 
8 The Nyquist critical frequency in the units of the plot is K, 20, 
well outside the range of v discussed in this section. 



corresponding to the second and third harmonics of V2 (i.e., 
2 V2 and 3 V2 , respectively) . The natural interpretation of 
this result based on the Fourier theorem is that /Edd is the 
composition of two periodic signals, with periodicities v\ and 

As an additional proof for the "physical" nature of the 
variability of /Edd, we computed an "a posteriori" estimate 
of /Edd by averaging the gas properties in the vicinity of the 
secondary every 2.5 x 10 4 yr and then the using these aver- 
ages to estimate the accretion rates. This approach allowed 
us to eliminate the effects of the sampling noise due to a 
finite number of SPH particles. More specifically, we com- 
puted the mean values of the density, relative velocity, and 
temperature of the gas inside the sphere of radius R p = 7 
pc centered around M2. The choice of the radius R p is mo- 
tivated by two arguments: R p is small enough to resolve the 
substructures in the disc caused by the orbiting MBH, and 
large enough to have a statistically large number of particles 
within R p from M%. The average number of particles inside 
the sphere defined by R p is ~ 20, 000, and the minimum 
number during the whole run is « 2, 000 (i.e., ~ 60 x -/V ne j g h). 
Such large numbers assure that the finite sampling of our 
disc does not affect our a posteriori estimate. Note that for 
counter-rotating M2, R P > -Rbhl- From these properties 
we estimated the Bondi-Hoyle-Littleton accretion rate onto 
M2. The power spectrum of the fluctuations of the new esti- 
mated /Edd is shown in Fig. 4 as a blue dot-dashed line. The 
positions of the second, third and fourth peaks are compati- 
ble with the V2 and its higher harmonics. The position of the 
first peak, different from the one found for our original data, 
can be explained considering that during the first phase of 
the first orbit M2 has not yet developed a parsec-scale over- 
density in the gas. Consequently, /Edd computed at large 
separations (i.e., within R p ) during the first orbit leads to 
lower signal in the frequency domain at v\ compared to the 
one obtained considering only bound particles (black curve). 
Note that despite this (expected) difference, the spectrum 
computed from the unaveraged data does possess a peak at 
v = v\, implying a physical origin of variability. 

In summary, the two tests discussed in this section prove 
the "physical" nature of the /Edd variability. The cause of 
such variability is discussed in the next section. 



4 DISCUSSION 

For the first time, we are able to study the dynamical evolu- 
tion of MBHs on retrograde orbits embeded in circumnuclear 
disks. These simulations require extremely high spatial res- 
olution (< 0.1 pc). We discovered that dynamical friction 
acting on a MBH moving on an initially retrograde orbit 
with respect to the disc material leads to an orbital angu- 
lar momentum flip. The moving MBH always ends up on 
a prograde orbit by the time the MBH pair forms a binary 
system (i.e., by the time the gas mass enclosed inside the 
MBH orbit is smaller than the sum of the MBH masses). 
We also find that the circularization of the initially eccen- 
tric prograde orbit of the secondary is efficient well before 
the formation of the binary. We stress that the orbital angu- 
lar momentum flip can only be captured provided that the 
numerical resolution is sufficiently high. 

The dynamical range of our simulations also allowed us 
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to resolve the sphere of influence of the MBHs and thus, for 
the first time, study the mass accretion rate of the MBH co- 
and counter-rotating with respect to the disc. We found that 
a highly variable double nuclear activity can be observed for 
few Myr when the two MBHs orbit each other with relative 
separations > 10 pc. The accretion history of the moving 
MBH (M2) can be divided into two distinct phases. The 
first stage occurs when the MBH is either on an eccentric 
co-rotating orbit or on a retrograde orbit irrespective of its 
eccentricity. The second stage occurs after the circulariza- 
tion of the orbit or after the orbital angular momentum flip 
of M2. In the first stage, the accretion rate of the secondary 
on a counter-rotating orbit is more variable, and its luminos- 
ity lower by a factor of a few, than in a co-rotating case. In 
the second stage, the accretion rate onto the moving MBH is 
much enchanced. Interestingly, the accretion rate of a MBH 
that is initially on a counter-rotating orbit exhibits a sig- 
nificantly more rapid transition to the high-mass accretion 
rate, low-variability state. The luminosity increases by an 
order of magnitude. 

We emphasise that the estimated accretion rate is com- 
puted on scales resolved by the simulations. We assume 
that the accretion rate is stationary. The validity of this 
assumption depends on the accretion disc physics on unre- 
solved scales, involving physical processes (i.e., magnetohy- 
drodynamics and microscopic transport processes) not im- 
plemented in our runs. A fully self-consistent study of the 
accretion onto the two MBHs down to few gravitational radii 
is beyond the scope of this investigation. 

The reason for the extreme variability in the earlier 
merger stage (out of the two described above) is twofold. 
Firstly, the mass accretion fluctuations are driven by the 
fast fluctuations in the relative velocity between M2 and 
the gas environment, that change the -Rbhl of M2, and by 
the inhomogeneities in the gas that passes in the vicinity of 
the secondary MBH. 

Secondly, in the counter-rotating cases, the finite num- 
ber of gas particles that we use to sample the gas distribu- 
tion may occasionally lead to short quiescence periods when 
the secondary is not accreting. In such situations, the -Rbhl 
of the secondary is not resolved for a brief time due to the 
extremly high relative velocity of the MBH with respect to 
the environment. The finite sampling implies also that val- 
ues of /Edd of order of 1CF 2 should be considered order of 
magnitude estimates, and can be incorrect by a factor of 
a few. However, values of /Edd larger than 0.1 have only a 
small relative numerical error of < 10 %. Therefore, the gen- 
eral dependencies of the luminosity variability on the initial 
merger parameters are robust. A MBH orbiting in a cir- 
cumnuclear disc during the final stages of a galaxy merger 
is expected to be embedded in a gaseous/stellar envelope 
comoving with the MBH itself. An accurate estimate of the 
amount of the gas and stars embedding the retrograde MBH 
depends on the initial conditions of the galaxy merger, and 
on the interplay between different physical processes, e.g., 
tidal and ram-pressure stripping. The presence of the stars 
and gas comoving with the MBH would increase the mass 
that perturbs the circumnuclear disc and, as a consequence, 
increase the effect of dynamical friction. This may reduce the 
timescale of the orbital angular momentum flip. An envelope 
of gas comoving with M2 can also increase the accretion rate 



onto the MBH, and decrease the variability of the accretion 
process, until the envelope is not totally accreted by M2 or 
removed by gravitational/ram-pressure stripping processes. 
We plan to study such effects in a future paper. 

The significant differences between /Edd for MBH on 
eccentric and circular orbits, and even stronger differences 
between / E dd in prograde and retrograde cases, can be con- 
sidered a possible explanation for the paucity of resolved 
double AGNs in circumnuclear discs. More specifically, it is 
conceivable that a circumnuclear disk in a galaxy that un- 
derwent a merger may host a pair of MBH, but only one of 
them may show clear signatures of an AGN. Alternatively, 
the secondary may be significantly dimmer, but highly vari- 
able. Such scenarios may occur when the M2 enters the disc 
on a counter-rotating trajectory. Such effects have impor- 
tant implications for predicting the rates of double AGN 
detections in galaxies. 

This study opens the way to predicting the properties of 
the accretion flows near the horizons of the coalescing MBH. 
Our simulations predict the amount of gas delivered to the 
very vicinity of the inspiraling MBH. As such our simula- 
tions serve as initial conditions required to model the proper- 
ties of the electromagnetic counterparts to the gravitational 
wave emission events that will be detectable with the Laser 
Interferometer Space Antenna (LISA) (Armitage & Natara- 
jan 2002; Kocsis et al. 2005; Milosavljevic & Phinney 2005, 
Dotti et al. 2006a; Lippai, Frei & Haiman 2008; Schnittman 
& Krolik 2008; Haiman, Kocsis & Menou 2008). We plan 
to improve upon our model by including gas cooling and 
heating, star formation and supernova feedback from newly 
formed stars. These new features can change the structure of 
the circumnuclear disc, creating a less homogeneous multi- 
phase environment (see, e.g., Wada & Norman 2001, Wada 
2004). The interaction between the two MBH and such a 
multiphase medium can in principle lead to a larger eccen- 
tricity of the forming binary, and to more variable accretion 
events onto the two MBH. 
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